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Abstract 



Pulsars are among the prime targets for the Large Area Telescope (LAT) aboard the 
recently launched Fermi observatory. The LAT will study the gamma-ray Universe 
between 20 MeV and 300 GeV with unprecedented detail. Increasing numbers of 
gamma-ray pulsars are being firmly identified, yet their emission mechanisms are 
far from being understood. To better investigate and exploit the LAT capabilities 
for pulsar science, a set of new detailed pulsar simulation tools have been developed 
within the LAT collaboration. The structure of the pulsar simulator package {Pul- 
sarSpectrum) is presented here. Starting from photon distributions in energy and 
phase obtained from theoretical calculations or phenomenological considerations, 
gamma rays are generated and their arrival times at the spacecraft are determined 
by taking into account effects such as barycentric effects and timing noise. Pulsars 
in binary systems also can be simulated given orbital parameters. We present how 
simulations can be used for generating a realistic set of gamma rays as observed by 
the LAT, focusing on some case studies that show the performance of the LAT for 
pulsar observations. 

Key words: Pulsars, gamma-ray pulsars, detectors, simulation 
PACS: 97.60.Gb, 95.85.Pw, 95.55.Ka 



1 Introduction 



Pulsars are among the most intriguing sources in the gamma-ray sky and 
represent unique laboratories for probing physical laws in extreme environ- 
ments. The EGRET experiment aboard the Compton Gamma Ray Obser- 
vatory (CGRO), whose mission ended in 2000, made the first all-sky survey 
above 50 MeV and made breakthrough observations of high-energy gamma- 
ray pulsars [39J. CGRO observed the already-known Vela and Crab pulsars, 
discovered five more gamma-ray pulsars, including the pulsed detection of 
Geminga, which at the time had been known just as a point source. Yet, 
many sources in the 3 rd EGRET catalog are unidentified [2fj, and pulsars are 
believed to account for the majority of unidentified sources near the Galactic 
plane [HlfTB] . The Fermi Gamma-ray Space Telescope, a new international 
and multi-agency space observatory, will be able to identify many of them 
unambiguously by searching for periodicity with high accuracy in space and 
time. 
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In spite of the efforts to understand the emission mechanisms in pulsars many 
questions still remain unanswered. A more detailed knowledge of high-energy 
gamma-ray emission from pulsars will help to unveil their real natures and the 
Fermi observatory, successfully launched aboard a Delta II rocket from Cape 
Canaveral Air Force Station on 11 June 2008, will provide deep insights into 
these fascinating objects. 

Fermi, formerly known as GLAST, carries two main instruments: the Large 
Area Telescope (LAT) which is a pair conversion telescope for gamma rays 
with energies between 20 MeV and 300 GeV, and the Gamma-ray Burst Mon- 
itor (GBM), devoted to the observation of Gamma Ray Bursts. 
The LAT consists of a high-precision silicon strip Tracker, a hodoscopic Csl 
Calorimeter, and a segmented Anticoicidence Detector which covers the full 
Tracker. A comprehensive description of the LAT can be found in [7j. 
Fermi will operate primarily in scanning mode, pointing away from the Earth 
and rocking about the zenith direction. The entire sky will be viewed every 
two orbits (~ 3 hours), providing long and uniform coverage. The high sensi- 
tivity due to the large effective area (~ 3800 cm 2 at E ~ 100 MeV and ~ 8000 
cm 2 for E>lGeV on-axis), the sharp Point Spread Function (# 68 ~0.1° for 
E>1 GeV), the large field of view (2.4 si) will allow the LAT to detect source 
fluxes down to ~ 3xl0~ 9 ph cm~ 2 s _1 This corresponds to a sensitivity 
~ 30 greater than EGRET's. An important factor for pulsar studies is the 
absolute timing accuracy of the LAT, which will be <10 /is for each detected 
gamma ray. Fermi has detected the EGRET pulsars and and so far detections 
of two new ones, the radio quiet pulsar in the CTA 1 supernova remnant [TJ 
and PSR J2021+3651 [2J, which has been also independently discovered by 
AGILE mission fl8l. 



To better understand the high-level performance of the LAT with regard to 
pulsar astronomy, a complete, new powerful simulation package (PulsarSpec- 
trum) has been developed. 

With respect to the simulation capabilities of existing packages, such as TEMP02 
[22J which have been used successfully for pulsar analysis by radio astronomers, 
Pulsar Spectrum is more focused on specific issues of interest in gamma-ray 
astronomy, as for example the possibility to include more complex spectra. 
Moreover, this simulator has the advantage of being specific for gamma-ray 
observations and fully compatible with the software that simulates the re- 
sponse of the LAT. 

The pulsar simulation package was a very useful tool before the launch of 
Fermi and it will continue to be so, for comparing Monte Carlo results with 
real data, making it easier to constrain the statistical significance of the re- 
sults. 

In pulsar simulations the most complex feature is to take into account any ef- 



for a steady source at E >100 MeV after 1 year sky survey and assuming an E 
spectrum with no spectral cutoff and a high latitude diffuse flux of 1.5 x 10~ 5 ph 
cm _2 s~ 1 sr~ 1 [7] 
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fects which can influence the photon arrival time at the telescope. They depend 
mainly on the spacecraft motion, and on the change of period due to rota- 
tional energy losses. Moreover, pulsars do not show a perfectly-steady increase 
of the period but exhibit timing noise, which becomes particularly important 
for young gamma-ray pulsars. This noise significantly reduces the phase co- 
herence with time and represents an important limitation when searching for 
pulsations over long time periods. PulsarSpectrum takes into account all these 
effects when computing the arrival times of the simulated gamma rays. It has 
been widely used by the LAT collaboration for full simulation of the expected 
pulsar population in the sky together with other sources. It has also been 
included in the LAT Standard Analysis Environment (SAE), that has been 
developed by the Fermi Science Support CenteiO} 



2 Gamma-ray pulsar simulation: motivations and basics 



In designing PulsarSpectrum all the observational information gathered so 
far about the known gamma-ray pulsars were used. EGRET discovered the 
gamma-ray pulsars PSR B1706-44 [38], PSR B1055-52 [13J and PSR B1951+32 
[12] by searching for gamma-ray pulsations at the same periods as their radio 
counterparts. Even though it was detected as a point source by SAS-2 and 
COS B, Geminga was not yet identified. The nature of Geminga as a gamma- 
ray pulsar was finally revealed when EGRET found a gamma-ray periodicity 
of about 237 ms starting from the X-ray pulsations detected by ROSAT [17] 
and identified Geming radio-quiet pulsaJ" 3 ! [8]. In addition, the BATSE 
and COMPTEL instruments aboard CGRO also detected PSR B1509-58 g3], 
even though this pulsar was never detected at EGRET energies. More recently, 
Fermi has discovered a new radio-quiet pulsar in the CTA 1 supernova rem- 
nant using blind search of periodicity in gamma rays pQ. 
The light curves of the EGRET detected pulsars are energy-dependent and 
above 100 MeV they show a double-peaked profile, while at > GeV energies 
one peak is reduced, usually the first peak coming after the radio pulse3- 
The peaks of the observed gamma-ray pulsars can be modeled as Lorentz pro- 
files (e.g. [23]). This profile has been chosen as a template when simulating 
random-generated light curves. 

Spectral analysis shows that most of the power output comes from the hard 
X-ray and gamma-ray emission, the latter being modeled by a non-thermal 
spectrum [ID]. EGRET spectral analysis of the brightest pulsars shows spec- 

2 http://fermi.gsfc.nasa.gov/ssc/ 

3 The Geminga pulsations were later found in COS B [9] and SAS-2 data [26] and 
the light curves were compared with the EGRET one[27j 

4 For all the simulations we adopt the convention of assigning phase to the radio 
pulse. 
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tral breaks at GeV energies, and observations by ACTs of the brightest pul- 
sars were unable to see pulsed emission [5j, although recently MAGIC has 
observed pulsed gamma rays above 25 GeV from the Crab pulsar [5]. The 
EGRET spectral observations and the absence of pulsed emission by ACTs 
(with the exception of the Crab) suggests a spectral cutoff that provides an 
important constraint on the main gamma-ray emission models proposed for 
pulsars [19]. In Polar Cap models, emission takes place at low altitude (< neu- 
tron star radius) above magnetic poles, and it is composed of emission from 
synchrotron-pair cascades that are initiated by curvature radiation converting 
via 7 — B —>■ e + e~ absorption [llj . Owing to the absorption in high B- fields, 
this model predicts sharp, super-exponential cutoffs in the observed spectra 
at ~ few GeV energies (See Eq. [[]). An evolution of the Polar Cap models 
are the Slot Gap models [29], where particle acceleration takes place in thin 
slot gaps along the last open field line from the neutron star surface to the 
light cylinder. In the Outer Gap model [33], instead, particles are accelerated 
within the vacuum outer gaps and emit gamma rays far from the neutron star 
surface mainly by curvature and synchrotron radiation. In high-altitude Slot 
Gap and Outer Gap models a simple exponential cutoff is predicted, due to 
the radiation reaction limit of the accelerated particles. Pulsar Spectrum has 
great flexibility for simulating different spectral shapes according to these the- 
oretical scenarios, since a power law spectrum with cutoff can be simulated 
in Pulsar Spectrum and its energy and shape can be adjusted according to the 
model. Thanks to this feature it will be possible to compare the distribution 
of expected photons with real data from the LAT, in order to constrain the 
theoretical emission scenarios. 

Phase- resolved spectral analysis of EGRET observations of Vela, Geminga and 
Crab pulsars show that no simple pattern can be easily recognized. However 
a trend in phase-dependent spectra has been found: the peaks appear softer 
than inter-peak emission and harder than the outer wings [31]. The depen- 
dence of the spectrum on phase is thus an important feature which needs to 
be simulated using a specific model in Pulsar Spectrum. 

In order to give some numerical estimates on the relative count rates obtained 
from simulations at different energy ranges, we used the Vela pulsar as an 
example, obtaining for E > 100 MeV a daily rate of ~ 440 counts from the 
pulsar, ~ 60 counts from the Galactic and extragalactic gamma-ray diffuse 
emission, ~ 5 counts from residual charged background and an upper limit of 
10 counts from the PWN [3]. At highest energies (E > 1 GeV), we obtained a 
daily rate of ~ 140 counts from the pulsar and ~ 16 counts from the Galactic 
gamma-ray diffuse, while the extragalactic diffuse and residual background 
give a negligible contribution. 
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3 Pulsar Spectrum overview 



Pulsar Spectrum generates a bi-dimensional gamma-ray photon histogram A^(E i; tj 
which represents the differential photon flux as a function of energy bins Ej_ 
and time bins U expressed in ph m -2 s -1 . N„ is used to extract energies and 
time intervals between the simulated photons, subsequently corrected for tim- 
ing effects (Sec. 13.31) . 

Pulsar Spectrum has two simulation strategies for creating N v : one is based on 
phenomenological considerations and uses analytical modeling of the pulsar 
spectrum; the other one is suited for simulations based on more complex pul- 
sar spectra, e.g. phase-dependent spectra that cannot be modeled analytically. 
After N v is created, PulsarSpectrum calculates the time interval between sub- 
sequent photons according to the source flux. Each generated photon is sent 
to the Monte Carlo simulation of the LAT for the reconstruction. 

3.1 The phenomenological model 

This model creates an N v flux histogram by multiplying a pulsar light curve 
and spectrum derived from phenomenology. The light curve can be randomly 
generated with one or two Lorentzian peaks or loaded from a file that lists the 
EGRET counts for each phase bin. Although a more detailed model should use 
the flux instead of counts, these are available only for the brightest EGRET 
pulsars [14] ; the light curves obtained are realistic enough for simulations. 
LAT light curves of faint EGRET pulsars are produced first by applying a 
smoothing boxcar filter to the EGRET light curves in order to reduce the 
statistical fluctuations, then by increasing the number of bins through linear 
interpolation to get time bin width of ~10 as, in accordance with the LAT 
requirements for absolute time resolution [35J. An example of this algorithm 
applied to the simulation of PSR B1055-52 is shown in Fig. [TJ Pulsar spectra 
are built by using the analytical formula of a power law with a super expo- 
nential cutoff [30]. This is a good model for the phase-averaged gamma-ray 
pulsar spectrum and allows for an adjustable shape of the energy cutoff. 
The parametric expression of the differential spectrum is: 



where a is the power law spectral index, Eq the cutoff energy, b the index 
which models the shape of the cutoff, and E n a scale factor set to 1 GeV for 
the standard simulations. These parameters are user selectable. 

5 In the following A^(Ej,tj) will be abbreviated with N u . 




bl 



(1) 
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Fig. 1. Example of creation of a simulated light curve based on EGRET observation 
of PSR B1055-52 using the phenomenological model in Pulsar Spectrum. The solid 
histogram reproduces the EGRET light curve from [40J. The dotted histogram is 
the result of a boxcar smoothing with 3-bins width. The solid curve is the smoothed 
histogram with increased number of bins and a phase shift of A0=O.l to assign the 
phase to the first bin. 

EqJH is useful for modeling the exponential cutoff (6=1) predicted by outer 
magnetosphere emission models or the super-exponential cutoff (b >1) pre- 
dicted by the Polar Cap emission model. The spectrum of EqJT] is multiplied 
by the light curve T(<f>) in order to obtain N u : 

where the rotational phases <fi have been converted to times t using the pulsar 
period PjjQ. The normalization factor C is computed in order to obtain an 
integrated flux above 100 MeV equal to the total flux of the pulsar in the 
EGRET energy band (100 MeV - 30 GeV) given in input. This model has been 
extensively used to check pulsar analysis tools and perform basic sensitivity 
studies [32J. 

3.2 The case of complex emission patterns 

For simulation of more complicated phase-energy distributions, another ap- 
proach must be used. The basic idea is simple: Pulsar Spectrum loads an ex- 
ternal N v bi-dimensional histogram and normalizes it according to the flux. 
In this way the computation of N u is delegated to external packages so many 
more scenarios can be used to produce a customized pulsar simulation. This 

6 For compatibility with the other source simulators developed by the LAT the 
time is used instead of the rotational phase. 
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kind of simulation is very useful to understand the utility of the LAT data for 
phase-resolved spectroscopy. As an example, we prepared a phase-dependent 
spectrum of the Vela pulsar based on EGRET observations [Hj . A super- 
exponential cutoff has been applied for the peak phase intervals and an expo- 
nential cutoff for the remaining phase intervals. From this model spectrum it 
is possible to estimate how the modeled light curve varies with energy (Fig. 
[2]). This model reproduces the behavior observed by EGRET, since the first 
peak is strongly reduced at high energies, leaving only the second one [4~T] . 
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Fig. 2. Light curves for the Vela pulsar at 3 energy ranges obtained from the 
phase-dependent spectral model based on EGRET observations. The spectral cutoff 
has been set to 8 GeV for the whole phase interval. The light curves for E>100 MeV 
and E>2 GeV have been magnified in order to compare the profiles in the same 
plot. 



3.3 Timing corrections 

The N u histogram is used to produce a list of photons, each one with assigned 
energy, arrival direction (from the pulsar position) and rotational phase. The 
simulator then calculates the arrival times as the LAT would measure them. 
In order to compute these arrival times, timing corrections must be applied 
so that after processing the simulated photons with the LAT pulsar tools 
the original light curve can be recovered. In the first stage the time intervals 
between photons are computed by assuming that: 

(1) the pulsar period is constant in time, 

(2) the LAT is fixed with respect to the pulsar, 

(3) there is no timing noise and 

(4) the pulsar is not in a binary system. 
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In this approximation the time intervals can be calculated considering the 
flux of the source and Poisson statistics. For fluxes comparable to those the 
gamma-ray pulsars, the time intervals computed by Pulsar Spectrum do not 
significantly differ from those computed for a pure Poisson process. The sim- 
ulator assume that the flux is constant with time, as observed from EGRET 
observations [23] . The assumptions listed above lead to simplified calculations 
but do not provide realistic time intervals St and so corrections must be ap- 
plied. In order to provide a realistic list of pulsar photons several effects must 
be incorporated as corrections to St, that are described below. 



3.3.1 Barycentric effects 

Since the Earth, and the Fermi observatory orbiting around it, undergo sig- 
nificant acceleration with respect to the pulsar, the first step in analyzing 
pulsar data is to convert the arrival times recorded at the spacecraft to the 
Solar System Barycenter (SSB) reference frame, which is nearly inertial with 
respect to the pulsar. 

Photon arrival times at the spacecraft are expressed in MET0Terrestrial Dy- 
namical Time (TT or TDT). The relation between the arrival time t#, ex- 
pressed in Barycentric Dynamical Time TDB0, and TT arrival time t F at the 
spacecraft is: 

t B = t F + At c {t F ) + At E (t F ) + At R (t F ) + At s (t F ) (3) 

where Atc(t F ) is the clock correction, At E (t F ) is the Einstein delay, Atn(t F ) 
is the Roemer delay and Ats{t F ) is the Shapiro delay. Starting from t%, the 
inverse of Eq. [3] is computed in order to give the arrival times t F at the 
spacecraft. 

Atc{t F ) is the correction from recorded time at the spacecraft to International 
Atomic Time and then to TT0. The Einstein delay At E (t F ) is a relativistic 
correction based on the variation in clock rate at the spacecraft and in the 
SSB because of the motion of the Earth and its gravitational potential. The 
values of Atc{t F ) and At F (t F ) are of ~1.5 /is and ~1.6 ms respectively [22] 



7 Mission Elapsed Time. For Fermi is defined as the number of seconds since 2001 
January 01 at 00:00:00 (UTC). 

8 According to IAU Standards the Barycentric Dynamical Time is the independent 
variable of the equations of motion with respect to the barycenter of the Solar 
System. It is related to TT by a mathematical expression that includes positions of 
Solar System bodies [22] . 

9 According to IAU Standards Terrestrial Time TT is the time reference for ap- 
parent geocentric ephemerides and it is related to the International Atomic Time 
(TAI) as as TT = TAI + 32.184 s [22]. 
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and are computed from the JPL ephemeridea I using the standard routine 
AXBARY.Q 
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The Roemer delay Atn(t F ) is due to the light propagation 
time from the position of Fermi to the SSB and it is based on the location 
of the Earth, the Sun and the Fermi spacecraft, and can reach ~500 s of 
amplitude. It can be written as: 



At R {t F ) = -- [r FE (t F ) + r ES (t F ) + T SB (t F )] ■ § 



(4) 



where r FE is the vector from Fermi spacecraft to the center of the Earth, 
r F s from the center of Earth to the center of the Sun and t$b from the center 
of the Sun to the SSB. s is the unit vector identifying the sky position of 
the pulsar, and can be considered constant, since proper motion and parallax 
effects are neglected by PulsarSpectrum.The Shapiro delay Ats{t F ) is a rela- 
tivistic correction due to the gravitational field of the Sun [31]. At the first 
order At s h(t F ) can be calculated as: 

2G 

M 8h (t F ) = log(l - s ■ f FS ) (5) 



where r F $ is the vector from Fermi to the Sun and G the gravitational con- 
stant. The Shapiro delay becomes significant only for sources that are located 
at small elongation from the Sun. Higher-order corrections are due to the Sun 
(~100/is), while corrections due to other planets (<180 ns) are neglected here 



3.3.2 Period changes 

Since the rotational energies of pulsars decrease with time, the periods in- 
crease. Pulsar Spectrum accounts for this effect by correcting the arrival times. 
When building a gamma-ray pulsar light curve the second step after barycen- 
tering is to assign a rotational phase to each photon. To do this, it is necessary 
to know the spin parameters, i.e. the rotational frequency f{to) and its deriva- 
tives /(to) an d /(to) a t a particular epoch t . In the inertial frame of the pulsar 
the rotational phase 0(t) can be written: 

(P(t) = frac[4>(t ) + /(t )(t - t ) + ^/(to)(t - t ) 2 + i/(t )(t - t ) 3 ]. (6) 

2 o 

where frac indicates the fractional part of the right hand member of the equa- 
tion. It is also possible to express the spin parameters as the period P and 

10 It is possible in Pulsar Spectrum to use the Jet Propulsion Laboratory DE200 or 
DE405 ephemerides |36j . 

11 See |http:/ /heasarc.gsfc.nasa.gov/docs/software/lheaso ft /RelNotes04.html 
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its derivatives P\ and P%. 

Higher-order terms in Eqj6j which better describe period evolution, are not 
believed to represent the secular spin down of the pulsar, but rather timing 
irregularities. For this reason they are included in the simulation of the timing 
noise. 



3. 3. 3 Timing Noise 

For most of the pulsars the evolution of the period does not follow a steady 
linear increase and fluctuations in phase can occur, because of the so-called 
timing noise. Since timing noise affects mainly young pulsars, it has a large 
impact on gamma-ray pulsars, which are mainly young objects. Thus a con- 
tinuous monitoring of radio and X-ray counterparts of gamma-ray pulsar can- 
didates is required [35]. Moreover, it reduces the phase coherence over time, 
and limits the search of periodicity using only gamma rays in Radio Quiet 
pulsars (blind search)^. 

The rotational phase <p(t) at an arrival time t can be modeled as a sum of 

a term 4>s(t), based on a steady increasing period (see Eq. E]), and a timing 

noise contribution </>tv(£). Pulsar Spectrum implements a routine for computing 

4>N(t) based on a Random Walk approach to timing noise |10j . 

Within this approach we model a Random Walk in the kth derivative of the 

phase: 

-^-± = J2a lU (t-U) (7) 



where u is the unit step function, tj are the times where discontinuous jumps 
in the fcth derivative occur, and are the amplitudes of the jumps. The 
interval between the times ti of two subsequent jumps is governed by a Poisson 
distribution with mean rate R. Within this Random Walk model, is a 
random variable with zero mean and Root Mean Square (RMS) that can be 
estimated from the pulsar spin parameters. According to the model developed 
in [TO] , Pulsar Spectrum is able to calculate the RMS in case of pure Phase 
Noise (k=0), pure Frequency Noise (k=l), pure Slowing- down Noise (k=2) 
or a linear superposition of them. 

Eventually, the RMS can be related to R, to the time T spanned by the data 
and to a term o"T'Ar(2,T), which represents the RMS residual phase from a 
least-squares 2 nd order polynomial fit over the interval T (See Eq. 21 in [10J). 
Summarizing, the relations used by Pulsar Spectrum to calculate the timing 
noise RMS are: 

r/2 °tnC2->T) 2 a^ N (2, T) 2 o~tn(2,T) 

dm oc — — ■ ; ov oc — — — ; ov oc — — — (8) 

Y RT RT 3 ' RT 5 V ) 
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In order to calculate the RMS for the different types of timing noise in Eq. 
IU Pulsar Spectrum uses T as the arrival time of the photon, while the mean 
rate R=l day -1 , in order to satisfy the validity condition RT >1 imposed by 
the model [TU]. The remaining term a TN (2, T) is computed using an activity 
parameter A, is defined as follows [10] : 

A = log w [a TN (m,T) / a TN (m,T) C rab]m=2 (9) 



where proportionality coefficients are taken from [TO], cr-r/vfm, T)c r qb repre 



sents the RMS of the phase residuals for the Crab pulsaa I and it is calcu- 



lated as in Eq. 16 of [10]. The first step is therefore to calculate the activity 
parameter, which is done by Pulsar Spectrum using the following relation [10] 



A = -1.37 + 0.71 log Pi 



(10) 



where the period derivative Pi is expressed in units of 10~ 15 s/s. Starting 
from the P of the simulated pulsar, Pulsar Spectrum calculates A (Eq. [TO]) . 
This value is used to evaluate o"Tjv(m, T) (Eq. [S]), and hence obtain 8(f) 2 , which 
represents the variance of the timing noise jumps (Eq. [H]). Fig. [^represents an 
example of simulated timing noise, which resembles the timing noise behavior 
shown in [TO]. This timing noise model is only one of the possible ones. An 



0.15 




150 200 

Time MET [days] 



Fig. 3. Example of the Random Walk timing noise implemented in PulsarSpectrum. 
In this case a pure phase timing noise has been applied to a simulated pulsar with 
period P = 0.28 s and a P\ = 1.1 x 10 -13 s/s, with a corresponding activity A~1.96. 

interesting possibility for a future implementation of PulsarSpectrum is to 



12 The RMS of the residuals for the Crab pulsar in Eq. [U]was chosen as the reference 
value because it was well-studied 1 10 1. 
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introduce the phenomenological timing noise model obtained from TEMP02 
timing solutions on specific pulsar observations. 



3.3-4 Orbital modulation and binary corrections 

Binary pulsars are an interesting target for the LAT. About 7% of the known 
pulsars 13 1 and almost 80% of millisecond pulsars (MSPs) are observed in bi- 
nary systems. Accretion from a companion star in fact is believed to be re- 
sponsible for the rejuvenation process that produces millisecond pulsars. 
In a binary system, the orbital motion and the gravitational field of the com- 
panion affect the timing of the pulsar and Pulsar Spectrum can simulate these 
effects starting from the orbital parameters. Dedicated routines for the correc- 
tion of these effects are included in LAT SAE. Orbits for non-relativistic binary 
systems can be described by Kepler's laws and calculations can be performed 
starting from 6 Keplerian parameters usually determined from observations 
at radio or other wavelengths [24] . 

The parameters required by Pulsar Spectrum to reference the photon arrival 
time to the barycenter of the binary system are the following: 

• orbital period 

• projected semi-major axis a p sini, where i is the orbital inclination, defined 
as the angle between the orbital plane and the plane of the sky; 

• orbital eccentricity e; 

• longitude of periastron u; 

• epoch of periastron passage T ; 

• position angle of the ascending node Q asc 

To compute the binary correction it is necessary to solve the Kepler equation, 
which can be written using the Eccentric Anomaly E(t) in the following way: 

E(t) -e sin E(t) = Q b [(t- T„)] (11) 



where = 2n/Pj, is the mean angular velocity. 

For every photon Pulsar Spectrum solves the Kepler equation numerically and 
finds the corresponding eccentric anomaly i?[23]. If the pulsar companion has 
a strong gravitational field, the description in terms of Keplerian parameters 
should be modified by including Post-Keplerian parameters. Since most of 
these parameters are very difficult to measure and are known for only a few 
pulsars |24| . only the Keplerian model is currently implemented in Pulsar- 
Spectrum. For pulsars in binary orbits the arrival time corrections described 

13 Source: ATNF Pulsar Catalog(http:/ / www.atnf.csiro.au/research/pulsar /psrcatj ) . 
Many simulations with Pulsar Spectrum made use of the ATNF Pulsar Catalog[25j. 
We would like to thank the ATNF team for maintaining this valuable resource. 
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in Eqj3] must be modified to include additional terms: 

t B = t F + At c + At E + At R + At s + At RB + At EB + At SB (12) 

Since General Relativistic effects are not fully implemented in the current 
version of Pulsar Spectrum At E B (the binary Einstein delay) and A$b (the 
binary Shapiro delay) are set to zero, and their magnitude is in the order of 
few [is [23]. The binary Roemer delay At R B across the pulsar orbit can be 
written in terms of the eccentric anomaly E as: 

Ausit) = x(cos E(t) — e) sintu + xsmE(t)yl — e 2 cosu; (13) 
where x = a p sini. 



4 Case studies of simulations 

In order to better understand the functionality of Pulsar Spectrum and its 
use in testing the analysis tools as well as the LAT capabilities for pulsar 
observations, some specific case studies are presented. 

4.1 PSR B 195 1+32 

The first example is a simulation of the light curve of PSR B1951+32, one 
of the faintest EGRET pulsars. The simulation was performed by assuming 
one year of observation in scanning mode. To model the gamma-ray spectrum 
of PSR B1951+32 we used a power law with spectral index a=1.74, a cutoff 
energy E =40 GeV and a super-exponential cutoff index 6=2 (spin param- 
eters have been taken from EGRET observations). The simulation included 
the diffuse gamma-ray background based on a model of cosmic rays and its 
interaction with the Interstellar Medium [37] developed by the LAT collabo- 
ration. For comparison with the previous EGRET results, photons within a 
1 degree circle around the source were selected. Simulated photons were then 
barycentered and phase-assigned using tools in the SAE. Following the same 
analysis steps performed on EGRET data, all the photons outside the phase 
intervals of the peaks (0.12< (p <0.22 and 0.48< <0.74) were assigned to 
the background and subtracted. The resulting light curve is shown in Fig. H] 
together with the one from EGRET. 

In this region EGRET collected 344 photons |42j. The statistics of the simu- 
lated LAT light curve are far superior to the EGRET observation, with the 
LAT light curve containing about 8 times as many gamma rays. Since the LAT 
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Fig. 4. Simulation of 1 year LAT observation of PSR B1951+32 in scanning mode 
compared with EGRET results. Top: 50 bin EGRET light curve for photons around 
1° from the pulsar position and energies above 100 MeV [22]. Bottom: 50 bin re- 
constructed light curve of LAT simulated data using the same region and energy 
range of EGRET data, compared with the simulated model (black line) . Both light 
curves are background-subtracted. 

has a wide energy range, a greater number of photons is also expected because 
of the hard spectrum of this source. This example suggests the potential of 
the LAT to discover new gamma-ray pulsars fainter than PSR B1951+32. 



4-2 Vela pulsar 



We also present study that illustrates the capability of the LAT to 

measure pulsar spectral cutoffs in order to distinguish between Polar Cap and 
Outer Gap models, as outlined in Sec. [TJ The Vela pulsar is the best candi- 
date for spectral cutoff measurements, since it is the brightest source in the 
gamma-ray sky and its cutoff energy should be around a few GeV, well within 
the LAT energy range. For this reason we simulated 1 year of LAT sky survey 
observations of the Vela pulsar using Polar Cap and Outer Gap models. The 
phase-averaged spectrum for Polar Cap is taken from the model of Daugherty 
and Harding 1996 [llj while the Outer Gap from the model of Romani 1996 
|33j . The diffuse gamma-ray background is the same as was used for PSR 
B1951+32 simulations. The data were analyzed with the XSpec vl2 pack- 



age I and a custom spectral model for fitting the super-exponential energy 



cutoff was used. Fig. [5] shows the resulting spectral distributions for the LAT 
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Fig. 5. Simulation of the Vela pulsar spectrum as observed in 1 year scanning mode. 
Polar Cap and Outer Gap spectrum are plotted together with the EGRET results 
to show the capability of the LAT to distinguish between them. 

with the EGRET data superimposed. From this plot is appears that LAT is 
sensitive enough to constrain the spectral model for Vela. 
It can be realistically estimated that the cutoff could be measured in a few 
months of observation in scanning mode, as it has been confirmed by Vela first 
observations by LAT{3]. However, to have a clearer insight into the details of 
the emission mechanism, a longer exposure is needed. The most powerful anal- 
ysis to study the emission mechanism is the phase-resolved spectral analysis, 
which should become feasible for the brightest pulsars within 1 year of data 
collection. 



4-3 Population Studies 



This example is based on the capability of simulating an entire pulsar pop- 
ulation, with the goal of studying the LAT sensitivity [32] or understanding 
the LAT detection ratio between radio-loud and radio-quiet pulsars [19J. We 
have produced large simulations based on a diversified and highly-detailed 
population built by using evolutionary synthesis codes. One example of such a 
population is shown in Fig. [6j where a total of 404 pulsars down to a flux limit 
of 10 _9 ph cm _2 s _1 (E>100 MeV) was produced starting from the pulsar dis- 
tribution provided by the population synthesis code described in [Tof2"0"] . The 
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Fig. 6. P-Pdot plot showing simulation of a population of 404 pulsars including 
different emission models, compared with the expected detections after 1 year (71 
normal pulsars and 61 MSPs). 

population is composed of the 6 EGRET pulsars and 39 radio pulsars coinci- 
dent with 3EG sources that were synthesized using the classical Polar Cap cas- 
cade scenario [15] . 140 pulsars with Low Altitude Slot Gap emission [28J were 
also simulated, by using a super exponential cutoff according to this model. 
103 were found to be fainter in the radio and were classified as radio-quiet, 
while the others were classified as radio- loud. The model in Fig. [6] includes also 
17 radio-loud and 212 radio-quiet millisecond pulsars (MSPs). According to 
the Polar Cap scenario, MSPs can produce gamma rays by Inverse-Compton 
induced cascades, even if they are below the Curvature Radiation pair death 
line [2H]. 

In order to estimate the number of pulsars detected with the LAT, we run 
an automatic routine based on the Fermi Science tools, obtaining a total of 
71 normal pulsars and 61 MSPs detected in 1 year. This result, that takes 
into account the detailed structure of the diffuse gamma-ray background, is 
compatible with the number of pulsars discovered discovered by Fermi in its 
first months [I]. 

The algorithms in Pulsar Spectrum have been optimized to keep simulation 
times short, and the simulation of large sample of pulsars have been used to 
estimate the computing time and precision. With a 2 GHz processor the sim- 
ulation of 200 pulsars for an observation time of 1 month takes nearly two 
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hours. The precision that PulsarSpectrum can achieve has a direct impact 
on the CPU time. In order to have a good compromise between precision and 
CPU time, the configuration adopted for PulsarSpectrum has an accuracy of 
50/is over the whole timing correction chain. This precision is considered sat- 
isfactory for the purpose of this simulation but can be improved if necessary. 



4-4 Timing noise and blind search 



The capability to simulate timing noise, which affects pulsar phase coherence 
over time (Sec. I3.3.3[) . allows a better understanding of the limits of the blind 
search for Radio Quiet pulsars. The first result of blind searches in the LAT 
data has been the discovery of the CTA 1 supernova remnant pQ. 
The timing noise produces a distortion in the light curve, as is shown in Fig. 
where an observation spanning 1 month is shown. With longer observation 
times the profile of the light curve continues to degrade until the pulsation is 
lost. The left panel shows a simulated pulsar located at CTA 1 with P =100 
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Fig. 7. Simulation of the effects of timing noise on a simulated pulsar over a month. 
Left: The pulsar is simulated without timing noise and phases are assigned using the 
exact values of Pq,P\, and Pi. Right: The same pulsar simulated including timing 
noise and with phase assignment using only an approximation of Poj and P\ as in 
a typical blind search. 

ms, Pi=10~ 13 s s _1 , and P2=9xl0~ 24 Hz ~ 2 and no timing noise. The phases 
have been assigned using the exact values of Po,Pi, and P 2 used for the simula- 
tion. The right panel shows a simulation of the same pulsar with timing noise, 
that has an activity parameter of A ~1.89 (Eq. [TO]) . The phase assignment 
has been carried out using approximate values of Pq, and Pi, as they would 
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be obtained from a blind search. In Fig. [7] a periodic structure is still visible 
but the significance is smaller, since the H-test statistics TS# change from 
TS H =209 (left) to TS H =125 (right). 

A more detailed set of simulations has been produced to explore how the 
timing noise RMS evolves with time for different activity parameters derived 
using typical values of P\ and keeping the same pulsar frequency of 10 Hz. 
We expect that higher activity parameters imply a more rapid increase of the 
timing noise RMS with time, as it is confirmed in Fig. [HJ An 0.1 ms or 1 ms 
RMS is compatible with typical RMS residuals provided from radio observa- 
tions. It is clear from Fig. [S]that a pulsar with a timing noise activity of 1.9 
reaches an RMS value of 0.1 ms over few days while the same value is reached 
in about 2 months for smaller timing noise. 
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Fig. 8. Evolution of the RMS of the simulated timing noise with time for different 
values of Pi corresponding to different values of activity A. 



5 Conclusions 



The first astronomical sources detected at gamma-ray energies were pulsars, 
but still many questions about their natures and emission mechanisms remain 
unanswered. Fermi provides an enormous leap in capability for exploring pul- 
sars and the LAT will be able to discover new gamma-ray pulsars, probing 
their nature with unprecedented detail. 

Simulation is a powerful tool to study the response of the instrument and 
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to check data analysis software. For this reason the LAT collaboration has 
produced a complete simulation package for the LAT detector, for different 
classes of sources and background, and a suite of software tools for analysis of 
gamma-ray data. 

Pulsar Spectrum, presented in this paper, is a new simulation package specific 
for gamma-ray pulsars. It has been demonstrated to be very flexible, accom- 
modating several alternative models, and very accurate, taking into account 
the main timing effects on the arrival time of the photons. The motion of the 
spacecraft in the Solar System with relativistic corrections and timing noise 
due to unpredictable changes in pulsar period and phase have been included. 
It has been demonstrated that timing noise introduces a substantial limitation 
when searching for radio-quiet pulsars. Pulsar Spectrum has been, and will re- 
main a very useful tool during the Fermi mission as a way to compare real 
data with the predictions of different theoretical models. 
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